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ABSTRACT 


We report the observation of a coalescing compact binary with component masses 2.5-4.5 Mo 
and 1.2-2.0 Mo (all measurements quoted at the 90% credible level). The gravitational-wave sig- 
nal GW230529_181500 was observed during the fourth observing run of the LIGO-Virgo-K AGRA 
detector network on 2023 May 29 by the LIGO Livingston observatory. The primary component of 
the source has a mass less than 5 Mo at 9996 credibility. We cannot definitively determine from 
gravitational-wave data alone whether either component of the source is a neutron star or a black 
hole. However, given existing estimates of the maximum neutron star mass, we find the most probable 
interpretation of the source to be the coalescence of a neutron star with a black hole that has a mass 
between the most massive neutron stars and the least massive black holes observed in the Galaxy. We 
estimate a merger rate density of 55r Gpc ? yr-! for compact binary coalescences with properties 
similar to the source of GW230529.181500; assuming that the source is a neutron star-black hole 
merger, GW230529_181500-like sources constitute about 60% of the total merger rate inferred for neu- 
tron star-black hole coalescences. The discovery of this system implies an increase in the expected rate 
of neutron star-black hole mergers with electromagnetic counterparts and provides further evidence 
for compact objects existing within the purported lower mass gap. 


Keywords: Gravitational wave astronomy(675) — Gravitational wave detectors(676) — Gravitational 
wave sources(677) — Stellar mass black holes(1611) — Neutron stars(1108) 


1. INTRODUCTION 


In May 2023, the fourth observing run (O4) of the Ad- 
vanced LIGO (Aasi et al. 2015), Advanced Virgo (Acer- 
nese et al. 2015), and KAGRA (Somiya 2012; Aso et al. 
2013) observatory network commenced following a series 
of upgrades to increase the sensitivity of the network. 
The prior three observing runs opened the field of ob- 
servational gravitational-wave (GW) astronomy, with 90 
probable compact binary coalescence (CBC) candidates 
reported by the LIGO Scientific, Virgo, and KAGRA 
Collaboration (LVK) at the conclusion of the third ob- 
serving run (O3) (Abbott et al. 2023a) and further can- 
didates found by external analyses (e.g., Nitz et al. 2023; 
Mehta et al. 2023b; Wadekar et al. 2023). These in- 
cluded the first observation of merging stellar-mass black 
holes (Abbott et al. 2016a), the first observation of a 
stellar-mass black hole merging with a neutron star (Ab- 
bott et al. 2021a), and the first observation of two merg- 
ing neutron stars (Abbott et al. 2017a). The first ob- 
servation of two merging neutron stars was also a mul- 
timessenger event accompanied by emission across the 


electromagnetic (EM) spectrum (Abbott et al. 2017b; 
Margutti & Chornock 2021). The continued discovery of 
CBCs by the Advanced GW network in O4 and beyond 
promises to reveal new information about the formation 
pathways of compact binaries and the physics of their 
evolution. 

While GW observations have enabled detailed charac- 
terization of the population of stellar-mass compact ob- 
ject binary mergers overall (e.g., Abbott et al. 2023b), 
the small number of observed mergers at the low-mass 
end of the black hole mass spectrum leads to consider- 
able uncertainty in this region of the mass distribution. 
Dynamical mass measurements of X-ray binary systems 
within the Milky Way suggest a paucity of compact ob- 
jects with masses between ~ 3-5 Mo, and hence a lower 
mass gap that divides the population of neutron stars 
(with masses less than ~ 3 Mo; Rhoades & Ruffini 1974; 
Kalogera & Baym 1996) and stellar-mass black holes 
(observed to have masses above 5 Mo; Bailyn et al. 
1998; Ozel et al. 2010; Farr et al. 2011b). Although a 
number of recent observations of non-interacting binary 
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systems (Thompson et al. 2019; Jayasinghe et al. 2021), 
radio pulsar surveys (Barr et al. 2024), and GW obser- 
vations (Abbott et al. 2020c, 2023b) have found hints 
of compact objects residing in the lower mass gap, GW 
observations have yet to quantify the extent and poten- 
tial occupation of this gap (Fishbach et al. 2020; Farah 
et al. 2022; Abbott et al. 2023b). 

Here we report on the compact binary merger 
signal GW230529_181500, henceforth abbreviated as 
GW230529, which was detected by the LIGO Livingston 
observatory on 2023 May 29 at 18:15:00 UTC; all other 
observatories were either offline or did not have the sen- 
sitivity required to observe this signal. We find that 
the compact binary source of GW230529 had component 
masses of 3.6* 05 Mo and 1.4*0$ Mo (all measurements 
are reported as symmetric 90% credible intervals around 
the median of the marginalized posterior distribution 
with default uniform priors, unless otherwise specified). 
Although we cannot definitively determine the nature 
of the higher-mass (primary) compact object in the bi- 
nary system, if we assume that all compact objects with 
masses below current constraints on the maximum neu- 
tron star mass are indeed neutron stars, the most prob- 
able interpretation for the source of GW230529 is the 
coalescence between a 2.5-4.5 Mo black hole and a neu- 
tron star. GW230529 provides further evidence that a 
population of compact objects exists with masses be- 
tween the heaviest neutron stars and lightest black holes 
observed in the Milky Way. Furthermore, if the source 
of GW230529 was a merger between a neutron star and 
a black hole, its masses are significantly more symmetric 
than the neutron star-black holes (NSBHs) previously 
observed via GWs (Abbott et al. 2023a), which increases 
the expected rate of NSBHs that may be accompanied 
by an EM counterpart. 

We report on the status of the detector network at the 
time of GW230529 in Section 2 and provide details of 
the detection in Section 3. In Section 4 we present es- 
timates of the source properties along with a discussion 
of the inferred masses, spins, tidal effects, and consis- 
tency of the signal with general relativity. We provide 
updated constraints on merger rates and the inferred 
properties of the compact binary and NSBH popula- 
tions as well as population-informed posteriors on source 
properties in Section 5. Section 6 provides analysis and 
interpretation of the physical nature of the source com- 
ponents. Implications for multimessenger astrophysics 
and the formation of low-mass black holes are in Sec- 
tions 7 and 8, respectively. Section 9 summarizes our 
findings. Data from the analyses in this work are avail- 
able on Zenodo (Abbott et al. 2024a). 


2. OBSERVATORY STATUS AND DATA QUALITY 


At the time of GW230529 (2023 May 29 18:15:00.7 
UTC), LIGO Livingston was in observing mode; LIGO 
Hanford was offline, having gone out of observing mode 
one hour and thirty minutes prior to the detection. The 
Virgo observatory was undergoing upgrades and was not 
operational at the time of the detection. The KAGRA 
observatory was in observing mode, but its sensitivity 
was insufficient to impact the analysis of GW230529. 
Hence, only the data from the LIGO Livingston obser- 
vatory are used in the analysis of GW230529. 

LIGO Livingston was observing with stable sensitiv- 
ity for ~ 66 hours up to and including the time of 
the GW signal. At the time of the detection, the sky- 
averaged binary neutron star (BNS) inspiral range was 
~ 150 Mpc. From the start of O4 until this event, the 
LIGO Livingston BNS inspiral range (Chen et al. 2021a) 
varied between 140-160 Mpc, a 4.5-19.4% increase com- 
pared to the median BNS range in O3 (Buikema et al. 
2020; Abbott et al. 2023a). Additional details on up- 
grades to the LIGO observatories for O4 can be found 
in Appendix A. 

The Advanced LIGO observatories are laser interfer- 
ometers that measure strain (Aasi et al. 2015). The ob- 
servatories are calibrated via photon radiation pressure 
actuation. An amplitude-modulated laser beam is di- 
rected onto the end test masses inducing a known change 
in the arm length from the equilibrium position (Karki 
et al. 2016; Viets et al. 2018). For the strain data used 
in the GW230529 analysis, the maximum 1e bounds on 
calibration uncertainties at LIGO Livingston were 6% 
in amplitude and 6.5° in phase for the frequency range 
20-2048 Hz. 

Detection vetting procedures similar to those of past 
GW candidates (Davis et al. 2021; Abbott et al. 2016b), 
when applied to GW230529, find no evidence that in- 
strumental or environmental artifacts (Helmling-Cornell 
et al. 2023; Nguyen et al. 2021; Effler et al. 2015) could 
have caused GW230529. We find no evidence of tran- 
sient noise that is likely to impact the recovery of the 
GW signal in the 256 s LIGO Livingston data segment 
containing the signal. 


3. DETECTION OF GW230529 


GW230529 was initially detected in low latency in 
data from the LIGO Livingston observatory. The sig- 
nal was detected independently by three matched-filter 
search pipelines: GSTLAL (Messick et al. 2017; Sachdev 
et al. 2019; Hanna et al. 2020; Cannon et al. 2021; Ew- 
ing et al. 2023; Tsukada et al. 2023), MBTA (Adams 
et al. 2016; Aubin et al. 2021) and PYCBC (Allen et al. 
2012; Allen 2005; Dal Canton et al. 2021; Usman et al. 


2016; Nitz et al. 2017; Davies et al. 2020). Although 
GW230529 occurred when only a single detector was 
observing, it was detected by all three pipelines with 
high significance, and stands out from the background 
distribution of noise triggers. More details are given in 
Appendix B. 

All three pipelines have a similar matched-filter-based 
approach to identify GW candidates but differ in the 
details of implementation. Each pipeline begins by per- 
forming a matched-filtering analysis on the data from 
the observatory with a bank of GW templates (Sakon 
et al. 2024; Roy et al. 2019; Florian et al. 2024). Times 
of high signal-to-noise ratio (S/N) are identified, after 
which each pipeline performs its own set of signal con- 
sistency tests, combines them with the S/N to make a 
single ranking statistic for each candidate GW event, 
and calculates a false alarm rate (FAR) by comparing 
the ranking statistic of the candidate with the back- 
ground distribution. The FAR is the expected rate of 
triggers caused by noise with a ranking statistic greater 
than or equal to that of the candidate. 

For each pipeline, this procedure can be done in two 
modes: a low-latency or online mode, and an offline 
mode. In the online mode, the pipelines only use the 
background information available at the time of a candi- 
date to estimate its significance, along with low-latency 
data-quality information. In the offline mode, pipelines 
use more background information gathered from sub- 
sequent observing times to estimate the significance of 
candidates and use more refined data-quality informa- 
tion. The offline mode enables more robust and repro- 
ducible results at the cost of greater latency. 

In both cases, the background distribution is extrap- 
olated to higher significances. This enables the inverse 
FAR to be greater than the time for which the back- 
ground was collected. Pipelines have differing extrapo- 
lation methods, and differing methods for calculating 
the FAR of a single-detector GW candidates, details 
of which can be found in Appendix C. These differ- 
ing methods can cause the FARs reported by different 
pipelines for the same candidate to vary to a signifi- 
cant degree, as seen in the case of GW230529. How- 
ever, all three pipelines recover a nearly identical S/N 
for GW230529, as expected for an astrophysical signal 
with a high match to the search templates. This, in con- 
cert with the fact that the inverse FARs from the offline 
analyses of all pipelines are much greater than the dura- 
tion of the first half of the fourth observing run (O4a), 
makes it highly likely that GW230529 is of astrophysi- 
cal origin. Table 1 shows the online S/N, online inverse 
FAR, and offline inverse FAR for each search pipeline. 
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Table 1. Properties of the detection of GW230529 from 
each search pipeline. Significance is measured by the inverse 
of the FAR. 


GstLAL MBTA  PyCBC 
Online S/N 11.3 11.4 11.6 
Online inverse FAR (yr) 1.1 1.1 160.4 


Offline inverse FAR (yr) 60.3 >1000 >1000 


4. SOURCE PROPERTIES 


The source properties of GW230529 are inferred us- 
ing a Bayesian analysis of the data from the LIGO Liv- 
ingston observatory. We analyze 128 s of data including 
frequencies in the range of 20-1792 Hz in the calculation 
of the likelihood. To describe the detector noise, we as- 
sume a noise power spectral density (PSD) given by the 
median estimate provided by BAYESWAVE (Littenberg 
& Cornish 2015; Cornish & Littenberg 2015; Cornish 
et al. 2021). We use the BiLBv (Ashton et al. 2019; 
Romero-Shaw et al. 2020) or PARALLEL BILBY (Smith 
et al. 2020) inference libraries to generate samples from 
the posterior distribution of the source parameters us- 
ing a nested sampling (Skilling 2006) algorithm, as im- 
plemented in the DvNESTY software package (Speagle 
2020). 

Given the uncertain nature of the compact objects of 
the source, we analyze GW230529 using a range of wave- 
form models that incorporate a number of key physi- 
cal effects. For our primary analysis, we employ binary 
black hole (BBH) waveform models that include higher- 
order multipole moments, the effects of spin-induced or- 
bital precession, and allow for spin magnitudes on both 
components up to the Kerr limit, but do not include 
tidal effects on either component. 

Systematic errors in inferred source properties due 
to waveform modeling may be significant for NSBH 
systems (Huang et al. 2021). We mitigate these ef- 
fects by combining the posteriors inferred using two dif- 
ferent signal models: the phenomenological frequency- 
domain model IMRPhenomXPHM (Pratten et al. 
2020; García-Quirós et al. 2020; Pratten et al. 2021), 
and a time-domain effective-one-body model SEOB- 
NRv5PHM (Khalil et al. 2023; Pompili et al. 2023; 
Ramos-Buades et al. 2023; van de Meent et al. 2023). 
'The posterior samples obtained independently using the 
two signal models are broadly consistent. Unless other- 
wise noted, we present results throughout this work ob- 
tained by an equal-weight combination of the posterior 
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Table 2. Source properties of GW230529 from the pri- 
mary combined analysis (BBH waveforms, high-spin, default 
priors). We report the median values together with the 
90% symmetric credible intervals at a reference frequency 
of 20 Hz. 


Primary mass mı /Mo 3.6108 
Secondary mass m2/Mo Prep 
Mass ratio q = m2/m1 Dagton 
Total mass M/Mo 5.1428 
Chirp mass M/Mo 1.9412:94 


Detector-frame chirp mass (1 + z)M/Mo 2.026+9:002 


. . ; 0.40 
Primary spin magnitude xi 0.444030 


Effective inspiral-spin parameter Yer —0.107 S47 
Effective precessing-spin parameter Xp DU 


02 
3015 


0.02 
0.04* 003 


Luminosity distance Dr,/Mpc 


Source redshift z 


samples from both models under the default priors de- 
scribed in Appendix D. Our measurements of key source 
parameters for GW230529 are presented in Table 2. 
Analysis details and results from other waveform mod- 
els we consider are reported in Appendix D; we find 
that the key conclusions of the analyses presented here 
are not sensitive to the choice of signal model. In par- 
ticular, the use of BBH models is validated by com- 
parison to waveform models that include tidal effects, 
finding no evidence that the BNS or NSBH models are 
preferred, consistent with previous observations (Abbott 
et al. 2021a). This is expected given the moderate S/N 
with which GW230529 was detected (Huang et al. 2021). 
The analysis of GW230529 indicates that it is an 
asymmetric compact binary with a mass ratio q — 
m»/m1 = 0.39+9:43 and source component masses mı = 
3.6193 Mo and ma = 1.4433 Mo. The primary is 
consistent with a black hole that resides in the lower 
mass gap (3 Mo S mi S 5 Mo; Ozel et al. 2010; 
Farr et al. 2011b), with a mass < 5 Mọ at the 99% 
credible level. The posterior distribution on the mass 
of the secondary is peaked around ~ 1.4 Mo with 
an extended tail beyond 2 Mọ, such that P(mz > 
2 Mo) = 5%. The mass of the secondary is con- 
sistent with the distribution of known neutron star 
masses, including Galactic pulsars (Antoniadis et al. 
2016; Alsing et al. 2018; Ozel & Freire 2016; Farrow 
et al. 2019) and extragalactic GW observations (Landry 
& Read 2021; Abbott et al. 2023b). Figure 1 shows 


GW170817 
GW190425 
GW190814 
GW?200105.162426 
GW200115_042309 


| GW230529 
0 10 20 
m; [Mo] 
Figure 1. The one- and two-dimensional posterior 


probability distributions for the component masses of the 
source binary of GW230529 (teal). The contours in the 
main panel denote the 90% credible regions with vertical 
and horizontal lines in the side panels denoting the 90% 
credible interval for the marginalized one-dimensional pos- 
terior distributions. Also shown are the two O3 NSBH 
events GW200105_162426 and GW200115_042309 (orange 
and blue respectively; Abbott et al. 2021a) with FAR < 
0.25 yr"! (Abbott et al. 2023a), the two confident BNS 
events GW170817 and GW190425 (pink and green respec- 
tively; Abbott et al. 2017a, 2019a, 2020a, 2024b), as well as 
GW190814 (red; Abbott et al. 2020c, 2024b) where the sec- 
ondary component may be a black hole or a neutron star. 
Lines of constant mass ratio are indicated by dotted gray 
lines. The grey shaded region marks the 3-5 Mo range 
of primary masses. The NSBH events and GW190814 use 
combined posterior samples assuming a high-spin prior anal- 
ogous to those presented in this work. The BNS events use 
high-spin IMRPhenomPv2 NRTidal (Dietrich et al. 2019a) 
samples. 


the component mass posteriors of GW230529 relative 
to other BNSs (GW170817 and GW190425) and NSBHs 
(GW200105.162426 and GW200115.042309, henceforth 
abbreviated as GW200105 and GW200115) observed by 
the LVK as well as GW190814 (Abbott et al. 2017a, 
2020a,c, 2021a). 

To capture dominant spin effects on the GW signal, 
we present constraints on the effective inspiral spin Veg, 
which is defined as a mass-weighted projection of the 
spins along the unit Newtonian orbital angular momen- 
tum vector Ly (Damour 2001; Racine 2008; Ajith et al. 


High spin secondary 
x1 < 0.99, x» < 0.99 


Low spin secondary 
T x1 « 0.99, x2 < 0.05 


Figure 2. Selected source properties of GW230529. 
The plot shows the one-dimensional (diagonal) and two- 
dimensional (off-diagonal) marginal posterior distributions 
for the primary mass mı, the mass ratio q, and the spin 
component parallel to the orbital angular momentum y1,z = 
Xic Ly. The shaded regions denote the posterior proba- 
bility with the solid (dashed) curves marking the 5096 and 
90% credible regions for the posteriors determined using a 
high spin (low spin) prior on the secondary of x2 < 0.99 
(x2 < 0.05). The vertical lines in the one-dimensional 
marginal posteriors mark the 90% credible intervals. 


2014), 


Xeft = (x T 2x2) Lu, (1) 
where the dimensionless spin vector x; of each compo- 
nent is related to the spin angular momentum S; by 
x; = cS;/(Gm2). If Xer is negative, it indicates that 
at least one of the spin component projections must be 
anti-aligned with respect to the orbital angular momen- 
tum, ie, Viz = Xi: Ly < 0. We measure an effective 
inspiral spin of Yer = —0.10*0-12, which is consistent 
with a binary in which one of the spin components is 
anti-aligned or a binary with negligible spins. The mea- 
surement is primarily driven by the spin component of 
the primary compact object y1,z = —0.11* 912. with a 
probability that y1,4 < 0 of 83%. However, there is a 
degeneracy between the measured masses and spins of 
the binary components such that more comparable mass 
ratios correlate to more negative values of Ye (Cutler & 
Flanagan 1994) for this system. We show the correlation 
between x1,, and the mass ratio and primary mass in 


Figure 2, with more negative values of x1,; correspond- 
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ing to more symmetric mass ratios and smaller primary 
masses. The secondary spin is only weakly constrained 
and broadly symmetric about 0, x5,; = —0.03*0 23. We 
find no evidence for precession, with the posteriors on 
the effective precessing spin yp (Schmidt et al. 2015) 
being uninformative. 

The presence of a neutron star in a compact binary im- 
prints tidal effects onto the emitted GW signal (Flana- 
gan & Hinderer 2008). The strength of this interaction is 
governed by the tidal deformability of the neutron star, 
which quantifies how easily the star will be deformed in 
the presence of an external tidal field. In contrast, the 
tidal deformability of a black hole is zero (Binnington 
& Poisson 2009; Damour & Nagar 2009; Chia 2021), of- 
fering a potential avenue for distinguishing between a 
black hole and a neutron star. We investigate the tidal 
constraints for both the primary and secondary com- 
ponents using waveform models that account for tidal 
effects (Dietrich et al. 2019a; Matas et al. 2020; Thomp- 
son et al. 2020a), which do not qualitatively change the 
mass and spin conclusions discussed above. Irrespec- 
tive of whether we analyze GW230529 with a NSBH 
model that assumes only the tidal deformability of the 
primary compact object to be zero or a BNS model that 
includes the tidal deformability of both objects, we find 
the tidal deformability of the secondary object to be 
unconstrained. The dimensionless tidal deformability of 
the primary peaks at zero, consistent with a black hole. 
The constraints on this parameter are also consistent 
with dense matter equation of state (EoS) predictions 
for neutron stars in this mass range. 

We also perform parametrized tests of the GW phase 
evolution to verify whether GW230529 is consistent with 
general relativity and find no evidence of inconsisten- 
cies. More detailed information on tidal deformability 
analyses and testing general relativity can be found in 
Appendix E.3 and Appendix F, respectively. 


5. IMPACT OF GW230529 ON MERGER RATES & 
POPULATIONS 


We provide a provisional update to the NSBH merger 
rate reported in our earlier studies (Abbott et al. 2021a, 
2023b) by incorporating data from the first two weeks 
of O4a using two different methods. In the first, event- 
based approach, we consider GW230529 to be represen- 
tative of a new class of CBCs and assume its contribu- 
tion to the total number of NSBH detections to be a sin- 
gle Poisson-distributed count (Kim et al. 2003; Abbott 
et al. 2021a) over the span of time from the beginning of 
the first observing run (O1) through the first two weeks 
of O4a. We find the rate of GW230529-like mergers to 
be R230529 = 55t Gpe? yr 1. When computing the 
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Figure 3. Posterior on the merger rates of NSBH systems. 
The solid and dashed lines represent the broad population- 
based rate calculation and the event-based rate calculation, 
respectively. 


rates of the significant NSBH events in O3 detected with 
FAR< 0.25 yr! (Abbott et al. 2021a, 2023a) using the 
same method, we find a total event-based NSBH merger 
rate of Rxspgu = 85*119 Gpc ? yr-!. Using this selec- 
tion criterion, the same as Abbott et al. (2023b), the 
population of NSBH mergers includes GW200105 and 
GW200115; we do not include GW190814, as its source 
binary is most probably a BBH (Abbott et al. 2020c; 
Essick & Landry 2020; Abbott et al. 2023b). 

For the second, broad population-based approach, we 
consider GW200105, GW200115, and GW230529 to be 
members of a single CBC class, together with an ensem- 
ble of less significant candidates. The definition of this 
class is determined by a simple cut on masses and spins 
resulting in a NSBH-like region of parameter space. We 
aggregate data from all the triggers found by GSTLAL 
from the beginning of O1 through the first two weeks of 
O4a, assess their impact on the estimated merger rates 
of NSBHs while also accounting for the possibility that 
some of them are of terrestrial origin (Farr et al. 2015; 
Kapadia et al. 2020), and update the NSBH merger 
rate obtained at the end of O3 (Abbott et al. 2023a) to 
RNSBH = 94+109 Gpc ? yr |. Further details regarding 
the classification of triggers in the population-based rate 
method can be found in Appendix G. 

We show updates to both the event-based and 
population-based NSBH rate constraints in Fig- 
ure 3. The event-based rate estimates highlight that 
GW230529-like systems merge at a similar (or poten- 
tially higher) rate to the more asymmetric NSBHs iden- 
tified in GWTC-3. The population-based rate estimate 
is more representative of the full NSBH merger rate as 


it includes less significant GW events in the data. Both 
of our updated estimates are consistent with the find- 
ings of Abbott et al. (2021a) within the measurement 
uncertainties. Further details about rate estimates can 
be found in Appendix G. 

In addition to updates to the overall merger rate of 
NSBHs, we also study the impact of GW230529 on the 
mass and spin distributions of the compact binary pop- 
ulation as inferred from GWTC-3 (Abbott et al. 2023b). 
We employ hierarchical Bayesian inference to marginal- 
ize over the properties of individual events and infer the 
parameters of a given population model (e.g., Thrane & 
Talbot 2019; Mandel et al. 2019; Vitale et al. 2020). The 
updates to population model results in this work are pro- 
visional because they only include one GW signal from 
O4a, although the biases resulting from this selection 
will not be severe since GW230529 occurred near the 
start of the observing run. We quantify this by compar- 
ing the event-based rate estimate (which accounts for 
O4a sensitivity to compute its detectable time-volume) 
with the NSBH rates attained by the various population 
analyses considered in this work, finding that they are 
consistent with each other. 

We use three different population models in our anal- 
ysis. The first two models consider the population of 
compact-object binaries as a whole without distinguish- 
ing by source classification, using either the parameter- 
ized POWER LAW + DiP + BREAK model (Fishbach 
et al. 2020; Farah et al. 2022; Abbott et al. 2023b) 
or the non-parametric BINNED GAUSSIAN PROCESS 
model (Mohite 2022; Ray et al. 2023b; Abbott et al. 
2023b). The POWER LAW + DIP + BREAK model is 
designed to search for a separation in masses between 
neutron stars and black holes by explicitly allowing for 
a dip in the component mass distribution. The BINNED 
GAUSSIAN PROCESS model is designed to capture the 
structure of the mass distribution with minimal assump- 
tions about the population. The broad CBC population 
analyses include all candidates reported in GWTC-3 
with FAR < 0.25 yr-!, the same selection criterion used 
in Abbott et al. 2023b (see Table 1 therein). The third 
model we investigate (NSBH-pop; Biscoveanu et al. 
2022) considers only the population of NSBH mergers 
with FAR < 0.25 yr™!. NSBH-pop is a parametric 
model designed to constrain the population distributions 
of NSBH masses and black hole spin magnitudes. This 
model assumes all analyzed events have a black hole 
primary and neutron star secondary; we do not include 
GW190814 in this analysis. Further details regarding 
the population model parameterizations and priors can 
be found in Appendix H. 
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Figure 4. The differential merger rate of NSBH systems as a function of black hole masses (solid curves: mean; shaded regions: 
90% credible interval) and minimum black hole mass (dashed lines: 90% credible interval) using the NSBH-POoP model with 
(teal) and without (grey) GW230529. The minimum black hole mass mpu,min is a parameter of the NSBH-POP population 
model, see Table 7. While the median rate is always inside the credible region, the mean can be outside the credible region. 


'The inclusion of GW230529 in our population analyses 
has several effects on the inferred properties of NSBH 
systems and the CBC population as a whole. 

The inferred minimum mass of black holes in 
the NSBH population decreases with the inclusion of 
GW230529. The mass spectrum of black holes in the 
NSBH population with and without GW230529 inferred 
using the NSBH-PoP model is shown in Figure 4. As 
the source of GW230529 is the NSBH with the small- 
est black hole mass observed to date, the minimum 
mass of black holes in NSBH mergers shows a sig- 
nificant decrease with the inclusion of this candidate: 
MminBH = 3.4773 Mo with GW230529 compared to 
Mmin,BH = 60155 Meg without. In contrast, the pa- 
rameter that governs the lower edge of the dip feature 
in the POWER LAW + DIP + BREAK model, which rep- 
resents the minimum black hole mass in the CBC pop- 
ulation, does not shift significantly with the inclusion 
of GW230529. This difference is because the POWER 
LAW + Dip + BREAK model makes no assumptions 
about the classification of the components and there- 
fore does not enforce sharp features at the edges of each 
sub-population, meaning that the source of GW230529 
is not necessarily a NSBH in this model. However, the 
BINNED GAUSSIAN PROCESS and POWER LAW + DIP 
+ BREAK models are designed to capture the structure 
of the full compact binary mass spectrum; because they 
do not assign a source classification to either of the bi- 
nary components, features present in these population 
models can have differing astrophysical interpretations 
from the NSBH-PoP model. 


GW230529 increases the inferred rate of compact bi- 
nary mergers with a component in the 3-5 Mo range. A 
region of interest in the mass distribution is the border 
between the masses of neutron stars and black holes. We 
choose ~ 3-5 Mo to represent the nominal gap between 
these populations. In Figure 5 we show the posterior on 
the rate of mergers with one or both component masses 
in the 3-5 Mo range, with and without GW230529. 
For the POWER LAW + Dip + BREAK model there 
is a small increase in the merger rate within this mass 
range, Rgap = 2475 Gpc ? yr-! with GW230529 ver- 
sus Rgap = 1717} Gpc ? yr! without. For the BINNED 
GAUSSIAN PROCESS model there is a larger increase 
in the merger rate, Rgap = 33153 Gpc ?yr-! with 
GW230529 versus Rgap = 7.5 554 Gpc ? yr! without. 
'The differing degree of change between the two models 
is due to different assumptions for the mass ratio distri- 
bution of merging compact binaries as well as the po- 
tential dip in the merger rate at low black hole masses 
built into the POWER LAW + DIP + BREAK model. 
BINNED GAUSSIAN PROCESS does not fit for a specific 
pairing function, while POWER LAW + Dip + BREAK 
assumes the same mass ratio distribution throughout 
the whole population of CBCs. As most observed BBHs 
favor equal masses, any population model conditioned 
on those observations should also favor equal masses in 
the BBH mass range. The implicit assumptions within 
the POWER LAW + DIP + BREAK model require this 
preference to be imposed for all masses, including rela- 
tively low-mass systems like GW230529. However, the 
assumptions in BINNED GAUSSIAN PROCESS do not 
broadcast this preference as strongly across different 
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Figure 5. Posterior on the merger rate of binaries with one 
or both components between 3-5 Mo. The solid curves show 
the results from the BINNED GAUSSIAN PROCESS analysis 
and the dashed curves show the results from the POWER 
LAW + Dip + BREAK analysis. Both models analyze the 
full black hole and neutron star mass distribution. The teal 
and grey curves show the analysis results with and without 
GW230529, respectively. 


mass scales and therefore may support more asymmet- 
ric mass ratios at lower masses. Regardless of the mass 
ratio distribution assumptions made by each model, we 
find that the inclusion of GW230529 provides further 
evidence that the ~ 3-5 Mo region is not completely 
empty. 

GW230529 is consistent with the population inferred 
from previously observed CBCs candidates. In Figure 6, 
we show the population distributions of the full black 
hole and neutron star mass spectrum for the primary 
component of compact binary mergers using the BINNED 
GAUSSIAN PROCESS (top panel) and POWER LAW + 
Dip + BREAK (bottom panel) population models. We 
qualitatively see that GW230529 is not an outlier with 
respect to the masses of previously-observed compact 
object binaries because the inclusion of GW230529 in 
the population does not significantly alter the full- 
population posterior constraints. This differs from the 
detection of GW190814, which was an outlier with re- 
spect to the rest of the observed BBH population at the 
time due to its small secondary mass (Essick et al. 2022). 
The observation of GW190814 strongly suggested the re- 
gion between the masses of neutron stars and black holes 
was populated (Abbott et al. 2023b), a conclusion that 
is strengthened with the detection of GW230529 (even 
though it is not an outlier). 

The component masses inferred for GW280529 differ 
across differing population-informed priors. In Figure 7, 
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Figure 6. The differential binary merger rate as a function 
of the mass of the primary component using the BINNED 
GAUSSIAN PROCESS model (top panel) and the POWER LAW 
+ Dip + BREAK model (bottom panel) for the full compact 
binary population (solid curves: mean; shaded regions: 90% 
credible interval) with (teal) and without (grey) GW230529. 


we show the mass and spin posteriors obtained using 
priors informed by each of the population models con- 
sidered in this work. The priors informed by the NSBH- 
POP model prefer unequal mass ratios and small black 
hole spins; they suppress the extended posterior tail out 
to equal masses and anti-aligned spins. The BINNED 
GAUSSIAN PROCESS model also pulls the posteriors to 
more asymmetric mass ratios and has less support for 
anti-aligned spins. As shown in the top panel of Fig- 
ure 6, the merger rate density inferred using the BINNED 
GAUSSIAN PROCESS model is nearly flat across the re- 
gion of parameter space covered by GW230529, mean- 
ing that the priors informed by this population model 
have less of an impact on the shape of the posterior 
compared to the parameterized models considered. Un- 
like the NSBH-POP model, the POWER LAW + DIP + 
BREAK model has a sharp drop in the merger rate above 
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Figure 7. Posterior distributions of GW230529 under vari- 
ous prior assumptions. The solid teal curve shows the poste- 
rior distributions using the default high-spin priors described 
in Appendix D. The various dashed and dotted curves show 
posteriors obtained using population priors based on the 
NSBH-pop model (orange), BINNED GAUSSIAN PROCESS 
(blue) and POWER LAW + Dip + BREAK (red) models. 


3 Mo. This sharp feature results in the posterior on the 
primary component being pulled below 3 Mo as well as 
the posterior on the secondary being pulled to higher 
masses. The BINNED GAUSSIAN PROCESS model has a 
similar feature, although it is closer to 2 Mo and there- 
fore too low to significantly affect the mass estimates for 
this system. The POWER LAW + DIP + BREAK model 
also assumes the same preference for equal-mass systems 
across the entire mass spectrum, and thus, given that 
the majority of the CBC population is consistent with 
symmetric mass ratios, infers that the GW230529 binary 
components are more similar in mass. The combination 
of the drop in the differential merger rate and the pref- 
erence for symmetric mass ratios increases the support 
for more extreme spins oriented in the hemisphere oppo- 
site the orbital angular momentum. From these models 
we find that the binary source of GW230529 either has 
asymmetric components with low values of veg, or has 
components that are similar in mass with xem ~ —0.3. 
'These results highlight that the choice of prior has a sig- 
nificant impact on the inferred masses and spins of the 
GW230529 source, and hence the inferred nature of the 
binary components. We assess the nature of the compo- 
nents further in Section 6. More details on the popula- 
tion prior reweighting are given in Appendix H.4. 


6. NATURE OF THE COMPACT OBJECTS IN 
THE GW230529 BINARY 


Without clear evidence for or against tidal effects in 
the signal, the physical nature of the compact objects in 
the GW230529 source binary can be assessed by com- 
paring the measured masses and spins of each compo- 
nent with the maximum masses and spins of neutron 
stars allowed by previous observational data. However, 
statistical uncertainties in component masses make it 
especially difficult to determine whether compact ob- 
jects with masses between ~ 2.5-5 Mo are consistent 
with being black holes or neutron stars (Hannam et al. 
2013; Littenberg et al. 2015). Nevertheless, we assess 
the nature of the source components by marginaliz- 
ing over our uncertainties in the masses and spins of 
GW230529, in the population of merging binaries, and 
in the supranuclear EoS, to compute the posterior prob- 
ability that each component had a mass and spin less 
than the maximum mass and spin supported by the 
EoS. We follow the procedure introduced in the context 
of GW190814 (Abbott et al. 2020c; Essick & Landry 
2020), which relied on only the maximum neutron star 
mass, and has subsequently been extended to include 
the effects of spin (Abbott et al. 2020c, 2021a, 2023b). 
Our analysis provides only an upper limit on the prob- 
ability that an object is a neutron star, as it assumes 
that all objects consistent with the maximum mass and 
spin of a neutron star are indeed neutron stars (Essick 
& Landry 2020). 

To assess the nature of the component masses of 
the GW230529 source, we consider two versions of an 
astrophysically-agnostic population model that is uni- 
form in source-frame component masses and spin mag- 
nitudes and isotropic in spin orientations: one with com- 
ponent spin magnitudes x; = |x;| that are allowed to be 
large (x1, X2 < 0.99) and one where they are restricted a 
priori to be small (x1, X2 € 0.05). We also consider the 
POWER LAW + Dip + BREAK population model (see 
Appendix H.2) fit to the GW candidates from GWTC- 
3 (Abbott et al. 2023a); including GW230529 in the 
population model does not significantly affect the re- 
sults. For each population, we marginalize over the un- 
certainty in the maximum neutron star mass conditioned 
on the existence of massive pulsars and previous GW ob- 
servations (Landry et al. 2020) using a flexible Gaussian 
process (GP) representation of the EoS (Landry & Es- 
sick 2019; Essick et al. 2020b). More information on the 
EoS choices can be found in Appendix I. 

'Table 3 reports the probability that each component 
of the merger is consistent with a neutron star. In gen- 
eral, we find that the secondary is almost certainly con- 
sistent with a neutron star, and the primary is most 
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Table 3. Probabilistic source classification based on consistency of component masses with the maximum neutron star mass 
and spin. All estimates marginalize over uncertainty in the masses, spins, and redshift of the source as well as uncertainty in the 
astrophysical population and the EoS. We consider three population models: two distributions that use astrophysically-agnostic 
priors and consider either large spins (x1, X2 < 0.99) or small spins (x1, x2 < 0.05), and a population prior using the POWER 
LAW + Dip + BREAK model fit with only the events from GWTC-3 (Abbott et al. 2023b). We use an EoS posterior conditioned 
on massive pulsars and GW observations (Landry et al. 2020). All errors approximate 90% uncertainty from the finite number 
of Monte Carlo samples used with the exception of the low-spin results for which we only place an upper or lower bound. 


X1, X2 < 0.99 X1,X2 < 0.05 POWER LAW + DIP + BREAK 
P(mı is NS) (2.9 + 0.4)% < 0.1% (8.8 + 2.8)96 
P(mz is NS) (96.1 + 0.4)% > 99.9% (98.4 + 1.3)% 


probably a black hole. However, when incorporating 
information from the POWER LAW + Dip + BREAK 
population model, we find that there is a ~ 1 in 10 
chance that the primary is consistent with a neutron 
star. If we further relax the fixed spin assumptions im- 
plicit within the POWER LAW + DIP + BREAK model 
for objects with masses € 2.5 Mo from x; < 0.4 to 
xi € 0.99 (see Appendix H.2), we can find probabilities 
as high as P(m, is NS) = (27.3 + 3.8)?6. This ambigu- 
ity is similar to the secondary component of GW190814 
(2.50 € m2/Mọo < 2.67), which is consistent with a neu- 
tron star if it was rapidly spinning (e.g., Essick & Landry 
2020; Abbott et al. 2020c; Most et al. 2020a). 

'The differences observed between population models 
primarily reflect the uncertainty in the mass ratio and 
spins of the GW230529 source. For example, incorporat- 
ing the POWER LAW + DIP + BREAK population model 
as a prior updates the posterior for m; from 3.6775 Mo 
to 2.773 Mo. Additional observations of compact ob- 
jects in or near the lower mass gap may clarify the com- 
position of GW230529 by further constraining the ex- 
act shape of the distribution of compact objects below 
~ 5 Mo or the supranuclear EoS. 


7. IMPLICATIONS FOR MULTIMESSENGER 
ASTROPHYSICS 


In NSBH mergers, the neutron star can either plunge 
directly into the black hole or be tidally disrupted by its 
gravitational field. Tidal disruption would leave some 
remnant baryonic material outside the black hole that 
could potentially power a range of EM counterparts in- 
cluding a kilonova (Lattimer & Schramm 1974; Li & 
Paczynski 1998; Tanaka & Hotokezaka 2013; Tanaka 
et al. 2014; Fernández et al. 2017; Kawaguchi et al. 2016) 
or a gamma-ray burst (Mochkovitch et al. 1993; Janka 
et al. 1999; Paschalidis et al. 2015; Shapiro 2017; Ruiz 
et al. 2018). The conditions for tidal disruption are de- 
termined by the mass ratio of the binary, the compo- 


nent of the black hole spin aligned with the orbital an- 
gular momentum, and the compactness of the neutron 
star (Pannarale et al. 2011; Foucart 2012; Foucart et al. 
2018; Krüger & Foucart 2020). While the disruption 
probability of the neutron star in GW230529 can be in- 
ferred based on the binary parameters, we are unlikely 
to directly observe the disruption in the GW signal with- 
out next-generation observatories (Clarke et al. 2023). 

We use the ensemble of fitting formulae collected in 
Biscoveanu et al. (2022) including the spin-dependent 
properties of neutron stars (Foucart et al. 2018; Cipol- 
letta et al. 2015; Breu & Rezzolla 2016; Most et al. 
2020b) to constrain the remnant baryon mass outside 
the final black hole following GW230529, assuming it 
was produced by a NSBH merger. We additionally 
marginalize over the uncertainty in the GP-EoS re- 
sults obtained using the method introduced in Sec- 
tion 6 (Legred et al. 2021, 2022). Using the high-spin 
combined posterior samples obtained with default pri- 
ors, we find a probability of neutron star tidal disruption 
of 0.1, corresponding to an upper limit on the remnant 
baryon mass produced in the merger of 0.052 Mo at 99% 
credibility. The low-secondary-spin priors (X2 < 0.05) 
yield a tidal disruption probability and remnant baryon 
mass upper limit of 0.042 and 0.011 Mo, respectively. 
A rapidly spinning neutron star is less compact than a 
slowly spinning neutron star of the same gravitational 
mass under the same EoS. This decrease in compactness 
leads to a larger disruption probability and a larger rem- 
nant baryon mass following the merger, explaining the 
trend we see when comparing the results obtained un- 
der the low-secondary- and high-spin priors. The source 
binary of GW230529 is the most probable of the con- 
fident NSBHs reported by the LVK to have undergone 
tidal disruption because of the increased symmetry in its 
component masses. However, the exact value of the tidal 
disruption probability and the remnant baryon mass for 
this system are prior-dependent. 
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Figure 8. Posterior on the fraction of NSBH systems de- 
tected with GWs that may be EM bright, fem-brignt, de- 
pending on the threshold remnant mass required to power 
a counterpart, J (Mes > ME cea) The solid and dashed 
curves represent different values of the minimum remnant 
mass Me us The teal and grey curves show the analysis 
results with and without GW230529, respectively. 


We can also gauge how the inclusion of GW230529 
impacts estimates of the fraction of NSBH systems de- 
tected in GWs that may be accompanied by an EM 
counterpart, fgw-brignt. Using the mass and spin dis- 
tributions inferred under the NSBH-POP population 
model described in Section 5, we find a 90% credi- 
ble upper limit on the fraction of NSBH mergers that 
may be EM-bright of fpm-bright < 0.18 when includ- 
ing GW230529, an increase relative to fEM-bright < 0.06 
obtained when excluding GW230529 from the analysis. 
'These estimates assume that any remnant baryon mass 
Mya min Z 0 could power a counterpart, although the 
actual threshold value is astrophysically uncertain. Fig- 
ure 8 shows the posterior for fEM-bright With and with- 
out the inclusion of GW230529. While the exact value 
of fem-bright depends on the assumed population model, 
the increase in fEM-bright for NSBHs upon the inclusion 
of GW230529 in the population is robust against mod- 
eling assumptions. 

Using these updated multimessenger prospects, we 
can infer the contribution of NSBH mergers to the pro- 
duction of heavy elements (Biscoveanu et al. 2022) and 
the generation of gamma-ray bursts (Biscoveanu et al. 
2023). Assuming that all remnant baryon mass pro- 
duced in NSBH mergers is enriched in heavy elements 
via r-process nucleosynthesis (Lattimer & Schramm 
1974, 1976), we infer that NSBH mergers contribute at 
most 1.1 Mc Gpc ? yr-! to the production of heavy 
elements and that the rate of gamma-ray bursts with 
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NSBH progenitors is at most 23 Gpc ? yr-! at 90% 
credibility. This likely represents a small fraction of 
all short gamma-ray bursts, for which astrophysical 
beaming-corrected rate estimates range from O(10 — 
1000) Gpc ? yr-! (Mandel & Broekgaarden 2022), and 
of the total r-process material produced by compact- 
object mergers (Coté et al. 2018; Chen et al. 2021b). 
No significant counterpart candidates have been re- 
ported for GW230529 (IceCube Collaboration 2023; 
Karambelkar et al. 2023; Lipunov et al. 2023; Longo 
et al. 2023; Lesage et al. 2023; Savchenko et al. 2023; 
Sugita et al. 2023a,b; Waratkar et al. 2023). This is un- 
surprising given that it was only observed by a single 
detector and hence was poorly localized on the sky. 


8. ASTROPHYSICAL IMPLICATIONS 


Since the late 1990s, there have been claims about the 
existence of a mass gap between the maximum neutron 
star mass and the minimum black hole mass (~ 3-5 Mo) 
based on dynamical mass measurements of Galactic X- 
ray binaries (Bailyn et al. 1998; Ozel et al. 2010; Farr 
et al. 2011b). The lower edge of this purported mass 
gap depends on the maximum possible mass with which 
a neutron star can form in a supernova explosion, which 
cannot exceed the maximum allowed neutron star mass 
given by the EoS. Some EoSs support masses up to 
~ 3 Mo for nonrotating neutron stars (Mueller & Serot 
1996; Godzieba et al. 2021) and even larger masses for 
rotating neutron stars (Friedman & Ipser 1987; Cook 
et al. 1994), although such large values are disfavored by 
the tidal deformability inferred for GW170817 (Abbott 
et al. 2019b) and Galactic observations (Alsing et al. 
2018; Farr & Chatziioannou 2020). The upper edge and 
extent of the mass gap depend on the minimum black 
hole mass that can form from stellar core collapse. How- 
ever, it remains an open question whether observational 
or evolutionary selection effects inherent to the detec- 
tion of Galactic X-ray binaries can lead to the observed 
gap in compact object masses (Fryer & Kalogera 2001; 
Kreidberg et al. 2012; Siegel et al. 2023). 

In recent years, new EM observations have unveiled 
a few candidates in the ~ 3 Mọ region, mostly from 
non-interacting binary systems (Thompson et al. 2019; 
van den Heuvel & Tauris 2020; Thompson et al. 2020b) 
and radio surveys for pulsar binary systems (Barr 
et al. 2024).  Microlensing surveys do not support 
the existence of a mass gap but cannot exclude it ei- 
ther (Wyrzykowski et al. 2016; Wyrzykowski & Man- 
del 2020). The LVK has already observed one compo- 
nent of a merger whose mass falls between the most 
massive neutron stars and least massive black holes 
observed in the Galaxy: the secondary component of 
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GW190814 (2.5-2.7 Mo at the 90% credible level; Ab- 
bott et al. 2020c, 2024b). The secondary components 
of the GW200210.092254 and GW190917_114630 source 
binaries also have support in the lower mass gap (Ab- 
bott et al. 2024b, 2023a), although their mass estimates 
are also consistent with a high-mass neutron star. Un- 
like GW230529, the primary components of all three 
of these binaries can be confidently identified as black 
holes. Overall, the existence of a mass gap between the 
most massive neutron stars and least massive black holes 
still stands as an open question in astrophysics. 

GW230529 is the first compact binary with a primary 
component that has a high probability of residing in the 
lower mass gap. Hence, GW230529 reinforces the con- 
clusion that the 3-5 Mo range is not completely empty 
(see Figure 5 in Section 5). However, the 3-5 Mo range 
may still be less populated than the surrounding regions 
of the mass spectrum. This conclusion is consistent with 
previous population analyses and rate estimates based 
on GWTC-3 (Abbott et al. 2023b). 

'The formation of GW230529 raises a number of ques- 
tions. Given our current understanding of core collapse 
in massive stars (O'Connor & Ott 2011; Janka 2012; 
Ertl et al. 2020), it is unlikely that the primary compo- 
nent formed via direct collapse because of its low mass, 
although stochasticity in the physical mechanisms that 
determine remnant mass may populate the low-mass end 
of the black hole mass spectrum (Mandel & Müller 2020; 
Antoniadis et al. 2022). Formation by fallback is a vi- 
able scenario: recent numerical models of core-collapse 
supernovae suggest that the formation of 3-6 Mọ black 
holes via substantial fallback is rare but still possible 
(e.g., Sukhbold et al. 2016; Ertl et al. 2020). One- 
dimensional hydrodynamical simulations of core collapse 
adopting pure helium star models predict that there is 
no empty gap, only a less populated region between 
3-5 Mo, with the lowest mass of black holes produced 
by prompt implosions starting at ~ 6 Mo (Ertl et al. 
2020). Another relevant parameter is the timescale for 
instability growth and launch of the core-collapse su- 
pernova; if this is long enough (= 200 ms), the proto- 
neutron star might accrete enough mass before the ex- 
plosion to become a mass-gap object (Fryer & Kalogera 
2001; Fryer et al. 2012). Compact binary population 
synthesis models that allow for longer instability growth 
timescales naturally produce merging NSBH systems 
with the primary component in the lower mass gap (Bel- 
czynski et al. 2012; Chattopadhyay et al. 2021; Zevin 
et al. 2020; Broekgaarden et al. 2021; Olejak et al. 2022). 
It may also be possible to form mass-gap objects through 
accretion onto a neutron star. If the first-born neu- 
tron star in the binary accretes enough material prior 


to the formation of the second-born compact object, it 
can trigger an accretion-induced collapse into a black 
hole, yielding a mass-gap object (Siegel et al. 2023). 
This scenario may be aided by super-Eddington accre- 
tion onto the first-born neutron star, which has been 
proposed to explain NSBHs with mass-gap black holes 
like the source of GW230529 (Zhu et al. 2023). Consid- 
ering the large number of uncertainties about the out- 
come of core-collapse supernovae (e.g., Burrows & Var- 
tanyan 2021), the primary mass of GW230529 provides a 
piece of evidence to inform and constrain future models. 
Overall, astrophysical models in the past have preferen- 
tially adopted prescriptions that enforce the presence of 
a lower-mass gap (e.g., Fryer et al. 2012), and the in- 
ferred rate of GW230529-like systems urges a change of 
paradigm in such model assumptions. 

Alternatively, the primary component of GW230529 
might be the result of the merger of two neutron stars. 
For instance, the 9096 credible intervals for the remnant 
mass of GW190425 and the primary mass of GW230529 
overlap (Abbott et al. 2020a). This may hint at a sce- 
nario where the primary component could either be the 
member of a former triple or quadruple system (Fra- 
gione et al. 2020; Lu et al. 2021; Vynatheya & Hamers 
2022; Gayathri et al. 2023), or the result of a dynami- 
cal capture in a star cluster (Clausen et al. 2013; Gupta 
et al. 2020; Rastello et al. 2020; Arca Sedda 2020, 2021) 
or AGN disk (Tagawa et al. 2021; Yang et al. 2020). 
'This scenario was proposed for the formation of a com- 
pact object discovered in a binary with a pulsar in 
the globular cluster NGC 1851, which is measured to 
have a mass of 2.09-2.71 Mo at 95% confidence (Barr 
et al. 2024). However, dynamically-induced BNS merg- 
ers are predicted to be rare in dense stellar environments 
(O(10-?) Gpc ? yr-!; e.g., Ye et al. 2020). Thus, the 
rate of mergers between a BNS merger remnant and a 
neutron star may be at least several orders of magnitude 
lower than the rate we infer for GW230529-like mergers, 
making this scenario improbable. 

Non-stellar-origin black hole formation scenarios such 
as primordial black holes (e.g., Clesse & Garcia-Bellido 
2022) remain a possibility. However, there are signifi- 
cant uncertainties in the predicted mass spectrum and 
merger rate of primordial black hole binaries, thus mak- 
ing it difficult to attribute a primordial origin to com- 
pact objects that have masses consistent with predic- 
tions from massive-star core collapse. Furthermore, re- 
sults from microlensing surveys indicate that primordial 
black hole mergers cannot be a dominant source of GWs 
in the local universe (Mroz et al. 2024). 

It has also been suggested that mergers apparently in- 
volving mass-gap objects could instead be gravitation- 


ally lensed BNSs (Bianconi et al. 2023; Magare et al. 
2023), with the lensing magnification making them ap- 
pear heavier and closer (Wang et al. 1996; Dai et al. 
2017; Hannuksela et al. 2019). This scenario is diffi- 
cult to explicitly test in the absence of tidal informa- 
tion (Pang et al. 2020) or EM counterparts (Bianconi 
et al. 2023), but the expected relative rate of strong lens- 
ing is low at current detector sensitivities (Smith et al. 
2023; Magare et al. 2023; Abbott et al. 2023c). 

Finally, we also find mild support for the possibility 
that the primary component is a neutron star rather 
than a black hole when considering the population-based 
POWER LAW + Dip + BREAK prior that incorporates 
the potential presence of a gap at low black hole masses. 
In this case, GW230529 would be the most massive neu- 
tron star binary yet observed, with both components 
= 2.0 Mo, and have non-zero spins that are anti-aligned 
with the orbital angular momentum. The effective in- 
spiral spin of the BNSs in this scenario would differ 
significantly from BNS sources previously observed in 
GWs (Abbott et al. 2017a, 2020a) as well as those in- 
ferred for Galactic BNSs if they were to merge (Zhu et al. 
2018). Spins oriented in the hemisphere opposite the bi- 
nary orbital angular momentum could be the result of 
supernova natal kicks (Kalogera 2000; Farr et al. 201 1a; 
O’Shaughnessy et al. 2017; Chan et al. 2020) or spin-axis 
tossing (Tauris 2022) at birth. For example, one of the 
neutron stars in the binary pulsar system J0737—3039 
is significantly misaligned with respect to both the spin 
of its companion and the orbital angular momentum of 
the binary system, which may require off-axis kicks (Farr 
et al. 2011a). Alternatively, isotropically-oriented com- 
ponent spins that result from random pairing in dynam- 
ical environments could lead to the significant spin mis- 
alignment (e.g., Rodriguez et al. 2016). 


9. SUMMARY 


GW230529 is a GW signal from the coalescence of a 
2.5-4.5 Mo compact object and a compact object con- 
sistent with neutron star masses. The more massive 
component in the merger provides evidence that com- 
pact objects in the hypothesized lower mass gap exist in 
merging binaries. Based on mass estimates of the two 
components in the merger and current constraints on 
the supranuclear EoS, we find the most probable inter- 
pretation of the GW230529 source to be a NSBH coales- 
cence. In this scenario, the source binary of GW230529 
is the most symmetric-mass NSBH merger yet observed, 
and the primary component of the merger is the lowest- 
mass primary black hole (BH) observed in GWs to 
date. Because NSBHs with more symmetric masses are 
more susceptible to tidal disruption, the observation of 
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GW230529 implies that more NSBHs than previously 
inferred may produce EM counterparts. However, we 
cannot rule out the contrasting scenario that the source 
of GW230529 consisted of two neutron stars rather than 
a neutron star and a black hole. In this case, the source 
of GW230529 would be the only BNS coalescence ob- 
served to have strong support for non-zero and anti- 
aligned spins, as well as the highest-mass BNSs system 
observed to date. Regardless of the true nature of the 
GW230529 source, it is a novel addition to the grow- 
ing population of CBCs observed via their GW emission 
and highlights the importance of continued exploration 
of the CBC parameter space in O4 and beyond. 

Strain data containing the signal from the LIGO Liv- 
ingston observatory is available from the Gravitational 
Wave Open Science Center.! Specifically, we release 
the L1: GDS-CALIB.STRAIN.CLEAN. AR channel, where the 
textttAR designation means the strain data is analy- 
sis ready; this strain channel was also released at the 
end of O3. Samples from the posterior distributions 
of the source parameters, hyperposterior distributions 
from population analyses, and notebooks for reproduc- 
ing all results and figures in this paper are available on 
Zenodo (Abbott et al. 2024a). The software packages 
used in our analyses are open-source. 
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Software: Calibration of the LIGO strain data 
was performed with GsTLAL-based calibration soft- 
ware pipeline (Viets et al. 2018). Data-quality prod- 
ucts and event-validation results were computed us- 
ing the DMT (Zweizig, J. 2006), DQR (LIGO Sci- 
entific Collaboration and Virgo Collaboration 2018), 
DQSEGDB (Fisher et al. 2021), gwdetchar (Urban et al. 
2021), hveto (Smith et al. 2011), iDQ (Essick et al. 
2020a), Omicron (Robinet et al. 2020) and PythonVir- 
goTools (Virgo Collaboration 2021) software packages 
and contributing software tools. Analyses in this cata- 
log relied upon software from the LVK Algorithm Li- 
brary Suite (LIGO Scientific, Virgo, and KAGRA Col- 
laboration 2018). The detection of the signals and sub- 
sequent significance evaluations were performed with 
the GsTLAL-based inspiral software pipeline (Messick 
et al. 2017; Sachdev et al. 2019; Hanna et al. 2020; 
Cannon et al. 2021), with the MBTA pipeline (Adams 
et al. 2016; Aubin et al. 2021), and with the Py- 
CBC (Usman et al. 2016; Nitz et al. 2017; Davies 
et al. 2020) packages. Estimates of the noise spectra and 
glitch models were obtained using BAYESWAVE (Cor- 
nish & Littenberg 2015; Littenberg et al. 2016; Cor- 
nish et al. 2021). Low-latency source localization was 
performed using BAYESTAR (Singer & Price 2016). 
Source-parameter estimation was primarily performed 
with the BILBy and PARALLEL BILBy libraries (Ashton 
et al. 2019; Smith et al. 2020; Romero-Shaw et al. 2020) 
using the DvNESTY nested sampling package (Spea- 
gle 2020). SEOBNRv5PHM waveforms used in param- 
eter estimation were generated using pySEOBNR, (Mi- 
haylov et al. 2023). PESUMMARY was used to postpro- 
cess and collate parameter-estimation results (Hoy & 
Raymond 2021). The various stages of the parameter- 
estimation analysis were managed with the Asimov li- 
brary (Williams et al. 2023). Plots were prepared with 
MATPLOTLIB (Hunter 2007), SEABORN (Waskom 2021) 
and GWPv (Macleod et al. 2021). NUMPY (Harris et al. 
2020) and SciPy (Virtanen et al. 2020) were used for 
analyses in the manuscript. 
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APPENDIX 


A. UPGRADES TO THE DETECTOR NETWORK FOR O4 


The Advanced LIGO, Advanced Virgo, and KAGRA detectors have all undergone a series of upgrades to improve 
the network sensitivity in preparation for O4. In this section, we focus on upgrades made to the LIGO Livingston 
observatory as it was the only detector to observe GW230529; similar upgrades were made at LIGO Hanford. Some 
of these improvements are a part of the Advanced LIGO Plus (A4-) detector upgrades (Abbott et al. 2020b). 

'The principal commissioning work done in preparation for O4 was to enhance the optics systems in the detector. 
The pre-stabilized laser was redesigned for O4 with amplifiers allowing input power up to 125 W (Bode et al. 2020; 
Cahillane & Mansell 2022). LIGO Livingston operated at 63 W in O4a. Both end test masses at LIGO Livingston 
were replaced because their mirror coatings contained small, point-like absorbers that limited detector sensitivity by 
degrading the power-recycling gain (Brooks et al. 2021). For O4, frequency-dependent squeezing was implemented 
with the addition of a 300 m filter cavity to rotate the squeezed vacuum quadrature across the bandwidth of the 
detector (Dwyer et al. 2022; McCuller et al. 2020). With frequency-dependent squeezing, uncertainty is reduced in 
both the amplitude and phase quadratures, allowing for a broadband reduction in quantum noise (Ganapathy et al. 
2023). Squeezing was implemented independently of frequency in O3, which only reduced the quantum noise at high 
frequencies (Tse et al. 2019). 

Other commissioning work was done to improve data quality and reduce noise at low frequencies. For O4, scattered 
light noise was mitigated by removing some problematic scattering surfaces, and improving the resonant damping of 
other scatterers (Soni et al. 2020; Davis et al. 2021). Low-frequency technical noise (i.e., noise that is not intrinsic to 
the detector design, but can limit the detector sensitivity and performance) was reduced with the commissioning of 
feedback control loops, noise subtraction, and better electronics. Overall, the upgrades made during the commissioning 
period for O4 led to a broadband reduction in noise and improvement in detector sensitivity and performance. 


B. DETECTION TIMELINE AND CIRCULARS 


The LVK issues low-latency alerts to facilitate prompt follow-up of GW candidates (Chaudhary et al. 2023). 
GW230529 was initially identified in a low-latency search using data from the LIGO Livingston observatory. The 
candidate was given the name $230529ay in the Gravitational-Wave Candidate Event Database (GRACEDB).? 

After the detection of GW230529, an initial notice was sent out to astronomers through NASA's General Coordinates 
Network (GCN) (LIGO Scientific, Virgo, and KAGRA Collaboration 2023a). The sky localization computed in low 
latency by BAYESTAR (Singer & Price 2016; Singer et al. 2016) had a 90% credible area of ~ 24200 deg?; the large 
credible area was due to GW230529 only being observed by a single detector. The initial alert also included low-latency 
estimates of the mass-based source classification (BNS, BBH, and NSBH) produced by the PYCBC pipeline (Villa- 
Ortega et al. 2022). Additional estimates that the source binary of GW230529 contained at least one neutron star, 
contained at least one compact object in the lower mass gap 3-5 Mo, and had neutron-star matter ejected outside the 
final compact object were produced by a machine learning method trained to infer source properties from the search 
pipeline results (Chatterjee et al. 2020). The estimated probability that the source binary of GW230529 included at 
least one neutron star was > 99.9%. 

Low-latency parameter estimation was performed with BILBY (Ashton et al. 2019; Romero-Shaw et al. 2020) and used 
to produce an updated sky localization and estimates of source properties sent out in another alert (LIGO Scientific, 
Virgo, and KAGRA Collaboration 2023b). The updated sky localization had a 90% credible area of ~ 25600 deg”. 
The updated estimate of source properties led to a slight decrease in the probability that the source binary included 
at least one neutron star, which was updated to 98.5%. 


C. QUANTIFYING SIGNIFICANCE FOR A SINGLE-DETECTOR CANDIDATE EVENT 


Each search pipeline has its own method for calculating the FAR of a single-detector GW candidate. In this section 
we will summarize these methods. 


? https: //gracedb.ligo.org/ 
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Figure 9. The noise background for the LIGO Livingston observatory collected by the GSTLAL pipeline in offline mode, with 
the S/N-£? values of GW230529 overlaid. The test statistic €? measures signal consistency. The background is collected from 
templates that are part of the same bin as the template that recovered GW230529, and is consequently the background used 
to rank GW230529 and calculate its significance. The background distribution has effectively no support at the position of 
GW230529. The colorbar represents the probability of producing a certain S/N and £? from noise triggers. 


C.1l. GsTLAL 


The GsTLAL pipeline assigns a likelihood ratio to all of its candidate signals as the ranking statistic to estimate 
their significance. The likelihood ratio is the ratio of the probability of obtaining the candidate under the signal 
hypothesis to the probability of obtaining the candidate under the noise hypothesis. The latter probability is largely 
calculated by collecting S/N and £? statistics of noise triggers during the analysis, where £? is a signal consistency 
test statistic (Messick et al. 2017). Similar templates in the template bank are grouped together on the basis of the 
post-Newtonian (PN) expansion of their waveforms (Morisaki & Raymond 2020; Sakon et al. 2024). Each template bin 
collects its own S/N and €? statistics, which get used to rank candidates recovered by templates in that bin. Further 
details about the calculation of the probability of obtaining the candidate under the noise hypothesis, as well as the 
GsTLAL likelihood ratio in general can be found in Tsukada et al. (2023). The S/N-£? statistics collected from noise 
triggers from the same template bin as GW230529 with GW230529 superimposed are shown in Figure 9. GW230529 
clearly stands out from the background, and there are no noise triggers at its position in S/N-£? space. 

Since non-stationary noise transients known as glitches (Nuttall 2018) are not expected to be correlated across 
detectors, single-detector GW candidates, whether during times when a single or multiple detectors are observing, 
are more likely to be glitches as compared to coincident GW candidates. To account for this, the GSTLAL pipeline 
downweights the logarithm of the likelihood ratio of single-detector candidates by subtracting an empirically-tuned 
factor of 13, which is optimized to maximize the recovery of candidates with an astrophysical origin while minimizing 
the recovery of glitches (Abbott et al. 2020a). The FAR of the GSTLAL pipeline for GW230529 quoted in Table 1 is 
calculated after the application of this penalty. 

Finally, the FAR for a candidate in the search is computed by comparing its likelihood ratio to the likelihood ratios 
of noise triggers not found in coincidence, after accounting for the live-time of the analysis. These noise triggers not 
found in coincidence allow the pipeline to extrapolate the background distribution of likelihood ratios to large values, 
enabling the inverse FAR of a candidate to be greater than the duration of the analysis. While in the low-latency 
mode, the GSTLAL FAR calculation uses the background from the start of the analysis period until the time of the 
candidate. In contrast, for its offline analysis, the GSTLAL pipeline uses background collected from the full analysis 
period including the entirety of O4a to rank candidates. This differs from the other search pipelines presented in 
this work which only use background from the first two weeks of O4a for their offline search results presented here. 
The GsTLAL offline analysis was performed using the same template bank as the online analysis (Sakon et al. 2024; 
Ewing et al. 2023). Details are liable to change for future offline GSTLAL analyses in O4a. However, for such a 
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Figure 10. Ranking statistic distribution of MBTA single-detector triggers in the LIGO Livingston observatory during the 
first two weeks of O4, excluding significant public alerts aside from GW230529. 


significant candidate as GW230529, we do not expect these changes to impact its interpretation as highly likely being 
of astrophysical origin. 


C.2. MBTA 


The MBTA single-detector candidate search for O4a focuses on a subset of the full MBTA parameter space. The goal 
is to focus on BNS and NSBH signals, which are most interesting because of their rarity and possible EM counterparts. 
Their long signal duration also makes them easier to identify and allows for a better sensitivity to be reached when 
compared to a single-detector candidate search using the whole MBTA parameter space. The parameter space for the 
O4a online single-detector search was defined based on the detector-frame (redshifted) primary mass and chirp mass 
of the templates and motivated by the computation of the probability for whether a non-zero amount of neutron star 
material remained outside the final remnant compact object (Chatterjee et al. 2020) introduced in Appendix B. The 
considered space is (1 + z) m, < 50 Mo and (1+ z) M <5 Mo (Juste 2023). For the offline analysis, this space was 
adapted and we consider only (1+ z) M < 7 Mo. MBTA online also applied selection criteria based on data-quality 
tests to its single-detector candidates (Juste 2023). This was changed to a reweighting of the ranking statistic for the 
offline analysis. 

The ranking statistic for MBTA single-detector triggers is a reweighted S/N. The reweighting is based on the 
computation of a quantity called auto x? (Aubin et al. 2021), which tests the consistency of the time evolution of the 
S/N time series. The additional reweighting used in the offline analysis relies on the computation of a quantity that 
identifies an excess of S/N for single-detector triggers. The excess of S/N is larger for triggers that have a noise origin 
compared to astrophysical signals or injections, and therefore allows for discrimination between astrophysical signals 
and noise. The ranking statistic distribution of the MBTA offline analysis for single-detector triggers during the first 
two weeks of O4a is shown in Figure 10. Other triggers which produced significant public alerts have been removed 
from the plotted distribution. GW230529 stands out from the background with a high ranking statistic that reflects 
the auto x? value being completely consistent with a signal origin. 

The FAR in MBTA is a function of Pastro (Andres et al. 2022), the probability of astrophysical origin of the 
candidate. It is derived from the combined parametrizations of the FAR as a function of the ranking statistic and 
Of Pastro as a function of the ranking statistic. Inverting the latter gives the ranking statistic as a function of Pastro 
and eventually the FAR as a function of Pastro. The background estimation for MBTA single-detector triggers was 
computed differently during the online and offline analyses. The online analysis relies on the computation of simulated 
single-detector triggers obtained through random combinations of single frequency-band data (Juste 2023). O4a online 
analysis and the method of randomly combining single-frequency band data showed that the background for MBTA 
single-detector triggers follows an exponential distribution and is stable in time beyond statistical fluctuations. This 
prompted us to update the model we use for the offline (and future) analyses to reach greater sensitivity. This new 
model involves extrapolating the observed distribution of single-detector triggers and removing the significant triggers. 
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Figure 11. Ranking statistic distribution of PyCBC offline single-detector triggers in the LIGO Livingston observatory during 
the first two weeks of O4. The plotted distribution may include triggers from other less-significant signals arriving during the 
period analyzed. 


A safety margin is used in the extrapolation such that the FAR is over-estimated relative to the best-fit extrapolated 
distribution. This change in methods, in addition to the difference in handling of single-detector triggers online and 
offline, explains the difference in inverse FAR for the online (1.1 yr) and offline (>1000 yr) analyses. 


C.3. PYCBC 


In PYCBC, each GW candidate is assigned a ranking statistic, and then a FAR is calculated by comparing the 
ranking statistic of the candidate to the background distribution. The details of this procedure differ between the 
online and offline versions of the PYCBC search. 

In the online pipeline, we consider single-detector candidates only from templates with duration greater than 7 s 
above a starting frequency of 17 Hz, corresponding to a range of masses and spins for which EM emission due to 
neutron star ejecta might be expected. Low-latency data quality timeseries produced by iDQ, a machine learning 
framework for autonomous detection of noise artifacts using only auxiliary data channels insensitive to GWs (Essick 
et al. 2020a), are used to veto candidates at times when a glitch is likely to be present in the data. The remaining 
single-detector candidates are ranked by the reweighted S/N described in Nitz (2018). Only candidates with a x? 
statistic (Allen 2005) < 2.0 and reweighted S/N > 6.75 are kept. The rate density of single-detector triggers above 
the reweighted S/N threshold is fit with a decreasing exponential. 

In offline PYCBC, we only consider single-detector candidates from templates with duration greater than 0.3 s 
above a starting frequency of 15 Hz. We also require that candidates have a x? statistic < 10.0, a reweighted S/N 
> 5.5, and a PSD variation statistic (Mozzon et al. 2020) < 10.0, in order to exclude high-amplitude noise transients. 
Each candidate is assigned a ranking statistic equal to the logarithm of the ratio of astrophysical signal likelihood to 
detector noise likelihood. The PyCBC O4 offline search introduced two changes to the ranking statistic relative to 
the O3 calculation. First, an explicit model of the signal distribution covering the space of binary masses and spins is 
included (Kumar & Dent 2024). 'The model is designed to maximize the number of detected signals from the known 
compact binary distribution, while also maintaining sensitivity to signals in previously unpopulated regions. Second, 
the model of the rate density of events caused by detector noise now includes a term describing the variation in rate 
during times of heightened detector noise (Davis et al. 2022). The rate density of single-detector candidates above a 
ranking statistic threshold of 0 is fit with a decreasing exponential. 

In both the online and offline versions of PYCBC, the exponential fit of trigger rate density above the relevant 
ranking statistic threshold is used to extrapolate the FAR of single-detector candidates beyond the observing time of 
the search (Cabourn Davies & Harry 2022). The FAR calculations for GW230529 used triggers from the first two weeks 
of O4 at the LIGO Livingston observatory. The distribution of ranking statistics for offline single-detector candidates 
at the LIGO Livingston observatory is shown in Figure 11. Similar to the other searches, GW230529 clearly stands 
out from the background distribution. 
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Table 4. Summary of parameter-estimation analysis choices for GW230529, and the physical content of each waveform model 
used. Here, tides refers to modeling of neutron star tidal deformability and disruption when tidal forces overcome the self-gravity 
of the neutron star. The spin prior denotes any restrictions on the spin magnitude of the binary components. 


Waveform Model Precession Higher Multipoles Tides Disruption Spin Prior 

IMRPhenomNSBH — — "4 "4 x1 < 0.50, x» < 0.05 
IMRPhenomPv2 NRTidalv2 "4 = Y — xi < 0.99, x2 < 0.05 
IMRPhenomXPHM "4 "4 — — xai < 0.99, x2 < 0.99 
SEOBNRv5PHM v "4 = = xai < 0.99, x2 < 0.99 
SEOBNRv4_-ROM_NRTidalv2_NSBH — — v v xai < 0.90, x2 < 0.05 
IMRPhenomXPHM "4 "4 — = xı < 0.99, x2 < 0.05 
IMRPhenomXP "4 = E — xi < 0.99, x» < 0.99 
IMRPhenomXHM — "4 — — xai < 0.99, x2 < 0.99 
IMRPhenomXAS — E = = xai < 0.99, x2 < 0.99 
IMRPhenomXAS — = = = xi < 0.50, x» < 0.05 
IMRPhenomPv2 NRTidalv2 "4 = v — xi < 0.05, x2 < 0.05 
IMRPhenomXPHM v — = xi < 0.05, x» < 0.05 
SEOBNRv5PHM "4 — — xai < 0.99, x2 < 0.05 


D. PRIORS, WAVEFORM SYSTEMATICS, AND BAYES FACTORS 


Given the uncertain nature of the compact objects, we analyze GW230529 with a suite of models that incorporate a 
number of key physical effects. The choices of data duration (128 s) and frequency bandwidth (20-1792 Hz) analyzed 
were informed by comparing waveforms spanning the mass range recovered in preliminary analyses to the detector 
noise PSD at the time of the event. The high-frequency cutoff is chosen to avoid loss of power at high frequencies due 
to low-pass filtering of the data. For a subset of the signal models below, we employ a range of techniques to speed 
up the evaluation of the likelihood including heterodyning (also known as relative binning; Cornish 2010; Zackay et al. 
2018; Cornish 2021; Krishna et al. 2023), multibanding (García-Quirós et al. 2021; Morisaki 2021), and reduced order 
quadratures (Canizares et al. 2015; Smith et al. 2016; Morisaki et al. 2023). 

'The physical effects included and the spin prior ranges for each waveform model considered in this work are shown 
in Table 4. All analyses use mass priors that are flat in the redshifted component masses with chirp masses, M = 
(mima)9/5 / (m, +m2)!/5 € [2.0214, 2.0331] Mo and mass ratios q € [0.125, 1]. The luminosity distance prior is uniform 
in comoving volume and source-frame time on the range Dy, € [1,500] Mpc. The priors on the tidal deformability 
parameters are chosen to be uniform in the component tidal deformabilities over A € [0, 5000]. Standard priors (e.g., 
Romero-Shaw et al. 2020; Abbott et al. 2024b) are used for all the other extrinsic binary parameters. Below we provide 
further motivation for considering this suite of waveform models. 

The primary analysis in this work uses the IMRPhenomXPHM and SEOBNRv5PHM BBH waveform models (corre- 
sponding to the third and fourth rows in Table 4), which provide an accurate description of BBHs but do not incorporate 
tidal effects or the potential tidal disruption of a companion neutron star. We did not consider an extension of IMR- 
PhenomXPHM that incorporates additional physics beyond the two precessing higher-order multipole moment BBH 
models used here as these effects only enter at high frequencies and are not relevant for this event (Thompson et al. 
2024). In order to quantify the impact of neglecting tidal physics, we analyze GW230529 with NSBH and BNS wave- 
form models that incorporate different tidal information. The NSBH models only incorporate tidal information from 
the less massive compact object, are restricted to spins aligned with the orbital angular momentum, and only model 
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the dominant £ = m = 2 multipole moment. However, they do model the possible tidal disruption of the neutron star, 
which occurs when tidal forces of the black hole dominate over the self-gravity of the neutron star. 

We use two NSBH models, a frequency-domain phenomenological model IMRPhenomNSBH (Thompson et al. 2020a), 
and a frequency-domain effective-one-body surrogate model SEOBNRv4 ROM NRTidalv2 NSBH (Matas et al. 2020). 
IMRPhenomNSBH uses the BBH models IMRPhenomC (Santamaria et al. 2010) and IMRPhenomD (Husa et al. 2016; 
Khan et al. 2016) as baselines for the amplitude and phase, respectively, incorporating corrections to the phase due 
to tidal effects following the NRTidal model (Dietrich et al. 2017, 2019b,a) and to the amplitude following Pannarale 
et al. (2013, 2015). SEOBNRv4_ROM_NRTidalv2_NSBH uses the SEOBNRv4 BBH model as a baseline (Bohé et al. 
2017) and applies the same corrections to account for tidal deformability and disruption as IMRPhenomNSBH but 
is additionally calibrated against numerical relativity (NR) simulations of NSBH mergers (Foucart et al. 2013, 2014, 
2019; Kyutoku et al. 2010, 2011). As the NSBH waveform models do not capture spin-induced orbital precession, 
we analyze GW230529 with the aligned-spin BBH waveform models IMRPhenomXAS (Pratten et al. 2020), which 
only contains the dominant harmonic, and IMRPhenomXHM (García-Quirós et al. 2020), which includes higher-order 
multipole moments. Models for NSBH binaries that contain both higher-order multipole moments and precession are 
under active development (Gonzalez et al. 2023) and could potentially allow for a more consistent model selection 
between the different source categories. 

Given that the mass of the primary overlaps with current estimates of the maximum known neutron star mass (Landry 
& Read 2021; Romani et al. 2022), we also analyze GW230529 with a BNS waveform model that allows for tidal interac- 
tions on both the primary and the secondary components of the binary. We use IMRPhenomPv2 NRTidalv2 (Dietrich 
et al. 2019a), which adds a model for the tidal phase that is calibrated against both effective-one-body (EOB) and NR 
simulations to the underlying precessing-spin point-particle waveform model IMRPhenomPv2 (Hannam et al. 2014). 
A limitation is that IMRPhenomPv2 NRTidalv2 is calibrated against a suite of equal-mass, non-spinning BNS NR 
simulations, where the maximum neutron star mass is only 1.372 Mo. Nevertheless, the model has been validated 
against a catalog of EOB-NR hybrid waveforms that includes asymmetric configurations and heavier neutron star 
masses (Dietrich et al. 2019a; Abac et al. 2024). Because of large uncertainties in BNS postmerger waveform modeling, 
IMRPhenomPv2 NRTidalv2 includes an amplitude taper from 1-1.2 times the estimated merger frequency (Dietrich 
et al. 2019b). The resulting suppression of S/N at these frequencies may be responsible for the posterior differences 
between IMRPhenomPv2 NRTidalv2 and the BBH and NSBH waveform models we consider (see Figure 12). 

From the mass estimates alone, the secondary object is consistent with being a neutron star. As such, we follow 
earlier analyses and analyze GW230529 with two spin priors (Abbott et al. 2021a): an agnostic high-spin prior and an 
astrophysically motivated low-spin prior. The high-spin prior assumes that the spins on both objects are isotropically 
oriented and have dimensionless spin magnitudes that are uniformly distributed up to x1, X2 < 0.99. The low-spin 
prior, on the other hand, is inspired by the extrapolated maximum spin observed in Galactic BNSs that will merge 
within a Hubble time (Burgay et al. 2003; Stovall et al. 2018), and restricts the spin magnitude of the secondary to 
be x2 < 0.05 while keeping xı < 0.99. As the nature of the primary as a black hole or neutron star is uncertain, we 
also use a third choice of prior where both spins are restricted to x1, x2 < 0.05 for the IMRPhenomPv2 NRTidalv2 
and IMRPhenomXPHM waveform models (rows 11 and 12 of Table 4, respectively). The purpose of using multiple 
priors is to help gauge whether astrophysical assumptions on the neutron star spin impact any of the statements on 
the probability that the primary object lies in the lower mass gap. In Figure 12, we show the impact of waveform 
systematics and prior assumptions on the strongly correlated mass ratio-effective inspiral spin distribution. 

The degree to which a given waveform model under certain assumptions matches the data can be gauged using 
the Bayes factor. However, Bayes factors penalize extra degrees of freedom in models that do not improve the fit 
when those extra degrees of freedom are constrained by the data (Mackay 2003). We find no significant preference 
(logio 8 X 0.5) between the high-spin and low-spin prior on the secondary; we similarly do not find a statistical 
preference for or against the effects of precession, higher-order multipole moments, or tidal deformation or disruption 
of the secondary object. 


E. ADDITIONAL SOURCE PROPERTIES 
E.1. Component Spins and Precession 


Assuming a high-spin prior, the posteriors for the primary spin magnitude are only weakly informative, disfavoring 
zero and extremal spins with xy, = 0.07-0.85. For the secondary, the posteriors are even less informative. Under 
the low-spin prior, the primary spin magnitude posterior is peaked at slightly higher values with zero and extremal 
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IMRPhenomXPHM 
xı < 0.99, x2 < 0.99 


IMRPhenomXPHM 
' xi < 0.99, x2 < 0.05 


SEOBNRv5PHM 
xı < 0.99, x2 < 0.99 


IMRPhenomNSBH 
xı < 0.5, v2 < 0.05 


SEOBNR_NSBH 
^ xi < 0.9, x» < 0.05 


IMRPhenomPv2_NRTidalv2 
x1 « 0.99, x» < 0.05 


Figure 12. The two-dimensional q—ye# posterior probability distributions for GW230529 using various BBH, NSBH, and BNS 
signal models. The vertical dotted lines indicate primary masses that have been mapped to the mass ratio given the median 
chirp mass estimated from the IMRPhenomXPHM high-spin analysis. 


spins being disfavored to a larger degree, xı = 0.10-0.84. The secondary spin is completely uninformative over the 
restricted range of spin magnitudes. The joint two-dimensional posterior probability distribution of the dimensionless 
spin magnitude and the spin tilt are shown in Figure 13. Regions of high (low) probability are denoted by a darker 
(lighter) shade. 

When the spins are misaligned with the orbital angular momentum, relativistic spin-orbit and spin-spin couplings 
drive the evolution of the orbital plane and the spins themselves (Apostolatos et al. 1994; Kidder 1995). The leading- 
order effect can be captured by an effective precession spin parameter, 0 < yp < 1, which approximately measures 
the degree of in-plane spin and can be used to parameterize the rate of precession of the orbital plane (Schmidt et al. 
2015) 


3+4 
Xp = max |a sin 61, (ica) q X2 sin A : (E1) 


We find the constraints on Xp to be uninformative, and we are not able to make any significant statements on 
precession. The uninformative nature of these results is corroborated by the Bayes factor between the precessing and 
non-precessing phenomenological waveform models, logio B2 = 0.22. 


E.2. Source Location and Distance 


As GW230529 was a single-detector observation, the sky localization is poor and covers a sky area of ~ 24100 deg? 
at the 90% credible level. The luminosity distance is inferred to be Dg = 201 Mpc, corresponding to a redshift 
of z= 0490.05 computed using the Planck 2015 cosmological parameters (Ade et al. 2016). The luminosity distance 
has a degeneracy with the inclination angle of the binary's total angular momentum with respect to the line of 
sight, Ozn (e.g., Cutler & Flanagan 1994; Nissanke et al. 2010; Aasi et al. 2013; Vitale & Chen 2018). The posterior 
distribution on zn is broadly unconstrained, showing no preference for a total angular momentum vector that is 


pointed towards or away from the line of sight. 
E.3. Tidal Deformability 
The tidal deformability of a neutron star is defined by 


2 
A= QR. (E2) 
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Figure 13. Two-dimensional posterior probability distributions for the spin magnitude x; and the spin-tilt angle 0; for the 
primary (left) and secondary (right) compact objects at a reference frequency of 20 Hz. A spin-tilt angle of 0? (180?) indicates 
a spin that is perfectly aligned (anti-aligned) with the orbital angular momentum L. The pixels have equal prior probability, 
and shading denotes the posterior probability of each pixel of the high-spin prior analysis, after marginalizing over azimuthal 
angles. The solid (dashed) contours denote the 9096 credible region for the high-spin (low-spin) prior analyses respectively. The 
probability distributions are marginalized over the azimuthal spin angles. 


where kg is the gravitational Love number of the object and R is its radius (Damour et al. 1992; Mora & Will 2004; 
Flanagan & Hinderer 2008; Damour & Nagar 2009). However, it is often convenient to introduce a dimensionless tidal 
deformability 


where C — m/R is the compactness of the neutron star in geometrized units. 

Using the IMRPhenomPv2 NRTidalv2 model, which allows for tidal deformability on both objects but does not 
model tidal disruption, we infer A; < 1462 at 9096 credibility with the x; < 0.99, x2 < 0.05 spin prior; the posterior 
peaks at A4 = 0. We do not constrain As relative to the prior. Constraints on the tidal deformability A of both 
components of GW230529 are shown in Figure 14, along with the tidal deformabilities predicted for the secondary 
object under the BSK24 neutron star EoS (Goriely et al. 2013; Pearson et al. 2018; Perot et al. 2019) and using 
the Gaussian process (GP)-EoS constraints introduced in Section 6. BSK24 is chosen as an illustrative EoS that is 
thermodynamically consistent and falls within the range of support of constraints from nuclear physics and astrophysics, 
including GW170817. The uninformative A» posterior is consistent with both of these predictions. In addition to the 
tidal deformability posteriors, we use two distinct likelihood interpolation schemes (Ray et al. 2023a; Landry & Essick 
2019) to directly constrain the neutron star EoS using both spectral (Lindblom 2010; Lindblom & Indik 2012, 2014) 
and GP (Landry & Essick 2019; Essick et al. 2020b) representations. Both methods incorporate consistency with the 
observed heavy pulsars PSR J0740--6620 (Fonseca et al. 2021) and PSR J0348+0432 (Antoniadis et al. 2013) as priors 
on the EoS and enforce thermodynamic stability and causality. Unlike other analyses where the GP representations 
are conditioned on both GW and pulsar data, here we only include the pulsar constraints in the prior, in order to 
distinguish the EoS information provided by GW230529 alone. These direct EoS inferences are also uninformative, 
returning their respective priors. 
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Figure 14. Constraints on dimensionless tidal deformability A for the primary (green dashed) and secondary (green solid) 
components of GW230529. We also show tidal deformabilities predicted for the secondary object under the BSK24 neutron 
star EoS (navy; Goriely et al. 2013; Pearson et al. 2018; Perot et al. 2019) and constraints obtained from BNS and pulsar 
observations using the GP model (GP-EoS; red). 


F. TESTING GENERAL RELATIVITY 


To verify whether GW230529 is consistent with general relativity, we perform parameterized tests searching for 
deviations in the PN coefficients that determine the phase evolution of the GW signal (Blanchet & Sathyaprakash 
1994, 1995; Arun et al. 2006a,b; Yunes & Pretorius 2009; Mishra et al. 2010; Li et al. 2012a,b). Specifically, we 
search for parametric deviations to the GW inspiral phasing applied to the frequency-domain waveform models IM- 
RPhenomXP_NRTidalv2 (Colleoni et al. 2023) and IMRPhenomXPHM (Pratten et al. 2020; García-Quirós et al. 
2020; Pratten et al. 2021) using the TIGER framework (Agathos et al. 2014; Meidam et al. 2018), and SEOB- 
NRv4_ROM_NRTidalv2_NSBH (Matas et al. 2020) and SEOBNRv4HM_ROM (Bohé et al. 2017; Cotesta et al. 2018, 
2020) using the FTI framework (Mehta et al. 2023a). These waveform models each capture different physical effects 
(precession, higher-order multipole moments, and neutron star tidal deformability) to determine whether their absence 
in the model leads to inferred inconsistencies with general relativity. 

For all waveform models and PN orders whose analyses have been completed, we find that GW230529 is consistent 
with general relativity within the inferred uncertainties on the deviation parameters. We do not yet include results 
about OPN deviations, which are being examined separately due to technical issues related to the degeneracy between 
the OPN deviation parameter and the chirp mass (Payne et al. 2023). The constraints obtained at —1PN are an order 
of magnitude tighter than previously reported bounds for NSBH and BBH (Abbott et al. 2021b). The previously 
reported bounds using GW170817 remain the tightest constraints on —1PN deviations obtained with GWs (Abbott 
et al. 2019c). Further details about the tests of general relativity performed will be presented in a follow-up paper. 


G. COMPACT BINARY MERGER RATES METHODS 


A key component of the event-based rate estimation described in Section 5 is the sensitive time-volume of our GW 
searches. The sensitivities to GW200105-, GW200115-, and GW230529-like events are computed from O1 through the 
first two weeks of O4a according to the corresponding mass and spin posteriors from the IMRPhenomXPHM high-spin 
analyses of these signals. The posterior samples chosen for this analysis are consistent with previous event-based rate 
estimates presented for NSBH mergers (Abbott et al. 2021a). Simulated signals, whose binary parameters are drawn 
from the mass and spin posteriors for each signal, are distributed uniformly in comoving volume and following the other 
extrinsic parameter distributions given in Appendix D, and are added to simulated detector noise characterized by 
representative PSDs for each observing run. The detectability of these injections is then calculated semi-analytically 
to determine the sensitive time-volume by calculating network responses to the simulated signals and applying a 
threshold on the network optimal S/N of > 10. This choice of threshold was previously tuned to the results of matched 
filter searches (Abbott et al. 2021c), and is comparable to threshold statistics calculated for semi-analytic sensitivity 
estimates (Essick 2023). 
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Table 5. Summary of BINNED GAUSSIAN PROCESS model parameters (Ray et al. 2023b; Mandel et al. 2017; Mohite 2022). 
Arguments in the priors specify the mean and standard deviation in a normal (N) or half-normal (HN) distribution. 


Parameter Description Prior 
Lu Mean log (a s) in each bin N(0, 10) 
o Amplitude of the covariance kernel HN(0, 10) 
log(l) log [e of the covariance kernel N(—0.085, 0.93) 


For the population-based approach, we estimate the merger rate of three astrophysical populations (BBH, BNS, and 
NSBH) by aggregating triggers found by GSTLAL from O1 through the first two weeks of O4, while accounting for 
the possibility that some of these triggers are of terrestrial origin. Our population analyses do not, however, include 
data from the engineering run preceding the start of O4, during which another NSBH candidate was identified (LIGO 
Scientific, Virgo, and KAGRA Collaboration 2023c). As outlined and implemented in previous LVK results (Abbott 
et al. 2021a, 2023b), we construct the joint likelihood on the Poisson parameters of the astrophysical populations 
(Annu, Apns, Anspu) and of the terrestrial triggers (Abackground; Farr et al. 2015; Kapadia et al. 2020). We then 
extrapolate the sensitive time-volume to each astrophysical population ((VT)ppu. (VT)nws, (VT) NspH) obtained at 
the end of O3 (Abbott et al. 2023b) to account for the additional two weeks of observation time during O4 during 
which GW230529 was detected. Using a uniform prior on the merger rates Ra = A4/(VT)s, we infer their posterior 
distributions. The three astrophysical populations are defined by dividing up the space of compact binary component 
masses and spins into disjoint regions. We consider all components with masses between 1 and 3 Mo to be neutron 
stars, and assume a maximum dimensionless spin of 0.05 for such components; all components that are more massive 
are assumed to be black holes with a maximum dimensionless spin of 0.99. The distribution of component masses 
within each region is assumed to be a power law in primary mass with index —2.35, and uniform in secondary mass. 

Even though we infer the posterior distributions of the merger rates of all three populations, we only present Rnspu 
given that GW230529 contributes negligibly to Rppy and Rpns. 


H. ADDITIONAL DETAILS OF POPULATION ANALYSES 
We use three different models to analyze the population of compact-object binaries with and without GW230529. 


H.1. BINNED GAUSSIAN PROCESS 


The BINNED GAUSSIAN PROCESS method models the merger rate density per log component masses as a piecewise 
constant function. By inferring the merger rate density in each bin, we reconstruct the shape of the CBC mass spectrum 
up to the resolution limit imposed by our choice of binning. A GP prior with an exponential quadratic kernel is imposed 
on the logarithmic rate densities to smooth out the inferred shapes over sparse regions of the parameter space. To 
assess the impact of GW230529 on the shape of the CBC mass spectrum, we use the same bin locations and priors on 
the GP hyper-paramters as the GWTC-3 analysis (Abbott et al. 2023b). The means (u), correlation length (l), and 
covariance amplitude (c) of the GP are drawn from normal, log-normal, and half-normal priors, respectively. Further 
details of these priors are summarized in Table 5. 

Unlike more recent implementations of the BINNED GAUSSIAN PROCESS model which also fit for the redshift dis- 
tribution (Ray et al. 2023b), we assume a redshift distribution such that the overall merger rate of compact binaries 
is uniform in comoving volume and source-frame time, so as to facilitate a direct comparison with the findings of the 
GW'TC-3 analysis (Abbott et al. 2023b) and to avoid any systematic biases in BINNED GAUSSIAN PROCESS-based 
redshift distribution models originating from the inclusion of low-mass events, which remains an active area of study. 
As in previous analyses (Abbott et al. 2023b; Ray et al. 2023b), we fix the spin distributions for each component to 
be isotropic in direction and uniform in spin magnitude. 
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Table 6. Summary of POWER LAw + Dip + BREAK model parameters (Farah et al. 2022). Arguments in the priors column 
specify the lower and upper bounds of a uniform (U) distribution. The first several entries describe the mass-distribution 
parameters, and the last two entries describe the spin-distribution parameters. 


Parameter Description Prior 
ay Spectral index for the power law of the mass distribution below My? U(—-8, 2) 
a2 Spectral index for the power law of the mass distribution above MAP U(—3, 2) 
A Lower mass gap depth U(0, 1) 
ME? [Mo] Location of the lower end of the mass gap U(1.4, 3) 
Msi; [Mc] Location of the upper end of the mass gap U(3.4, 9) 
Mow Parameter controlling how the rate tapers at the low end of the mass gap 50 
T]high Parameter controlling how the rate tapers at the high end of the mass gap 50 
n Parameter controlling tapering the power law at high black hole mass U(—4, 12) 
B Spectral index for the power law-in-mass-ratio pairing function U(-2, 7) 
Mmin [Mo] Minimum mass of the mass distribution U(1, 1.4) 
Mmax [Mc] Maximum mass of the mass distribution U(35, 100) 
Xmax,NS Maximum allowed component spin for objects with mass < 2.5Mo 0.4 
Xmax,BH Maximum allowed component spin for objects with mass > 2.5Mo 0.99 


H.2. POWER LAW + Dip + BREAK 


The POWER LAW + Dip + BREAK model is designed to search for a separation in masses between neutron stars and 
black holes by employing a broken power law with a dip at the location of the power-law break. The dip is modeled 
by a notch filter with depth A, which is fit along with other model parameters in order to determine the existence and 
depth of a potential mass gap (Farah et al. 2022). A value A = 0 corresponds to no gap, whereas A = 1 corresponds to 
zero merger rate over the interval of the gap, i.e., a maximally deep gap. POWER LAW + DIP + BREAK also employs 
a low-pass filter at high black hole masses to allow for a tapering of the mass spectrum, which has the effect of adding 
a smooth second break to the power law. 

The component mass distributions in this model are both fit by the same broken power law with exponents o4 
between Mmin and MP and az between MEY and mmax- The model additionally includes a power-law pairing 
function in mass ratio (Fishbach & Holz 2020), assumed to be the same for all component masses. The parameters for 
this mass model are summarized in Table 6. Like the BINNED GAUSSIAN PROCESS model, the POWER LAW + DIP + 
BREAK model additionally assumes that CBCs are uniformly distributed in comoving volume and source-frame time 
and that component spins are isotropically oriented and uniformly distributed in magnitude. These assumptions are 
made for simplicity and consistency with GWTC-3 analyses (Abbott et al. 2023b). Components with m < 2.5 Mo 
are limited to spin magnitudes < 0.4 while components with m > 2.5 Mọ can have spin magnitudes up to 0.99. 


H.3. NSBH-PoP 


The NSBH-PoP model introduced in Biscoveanu et al. (2022) is designed to capture the mass and spin distributions 
of the NSBH population. The black hole mass distribution is fit by a truncated power law, and the conditional mass 
ratio distribution, p(q|mi), is fit by a truncated Gaussian between qmin = 1 Mo /m and qmax = min(mws,max/m, 1). 
The black hole spin magnitude distribution is modeled as a Beta distribution (including singular distributions that 
peak at the edges of the x1 € [0,0.99] space), while the neutron star spin magnitude distribution is assumed to be 
uniform over x» € [0, 0.7] to encapsulate the effect of the neutron star breakup spin at the mass-shedding limit (Shao 
et al. 2020; Most et al. 2020a). We assume the redshift distribution is uniform in comoving volume and source-frame 
time and that spins are isotropically oriented. The parameters for this model are summarized in Table 7. 
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Table 7. Summary of NSBH-pop model parameters (Biscoveanu et al. 2022). Arguments in the priors column specify the 
lower and upper bounds of a uniform (U) distribution. 


Parameter Description Prior 
a Black hole mass power-law index U(—4, 12) 
MBH,min [Mc] Minimum black hole mass U(2, 10) 
MBH,max [Mo] Maximum black hole mass U(8, 20) 
MNS,max [Mo] Maximum neutron star mass U(1.97, 2.7) 
m Mass ratio mean U(0.1, 0.6) 
c Mass ratio standard deviation U(0.1, 1) 
Ox Beta distribution shape parameter (a) for black hole spin distribution U(0.1, 10) 
Bx Beta distribution shape parameter (3) for black hole spin distribution U(0.1, 10) 


H.4. Population reweighting 


We use the population-level mass and spin distributions inferred using all three population models to reweight 
(e.g., Payne et al. 2019) the high-spin combined parameter estimation samples from the default priors described in 
Appendix D to three different astrophysically-motivated priors given by the one-left-out posterior predictive distribu- 
tions (Galaudage et al. 2020; Callister 2021; Essick & Fishbach 2021) inferred for each model. A comparison of the 
posteriors under the default and astrophysically-motivated priors is shown in Figure 7. 

Because of the q— Xer degeneracy, the spin distribution assumptions also impact the inferred component masses, and 
hence classification, for this source. However, the fact that the BINNED GAUSSIAN PROCESS and NSBH-PoP model 
priors recover similar mass posteriors indicates that the different spin distribution assumptions have a subdominant 
effect compared to the different mass distribution assumptions. 


H.5. Selection effects 


For our population analyses, we account for search selection bias to measure the underlying astrophysical mass and 
spin distributions of the NSBH and CBC populations (Loredo 2004; Mandel et al. 2019; Farr 2019; Vitale et al. 2020). 
To estimate the sensitivity of the searches over the data recorded by the detector network across the binary parameter 
space, we use a large suite of injections and impose a significance threshold on the FAR. In contrast to the two 
methods for calculating merger rates described in Appendix G, we use the same set of injections used for GWTC-3 
population studies recovered using matched-filter searches (Abbott et al. 2023b). By using the combined injection sets 
from the first three observing runs (Abbott et al. 2021c), we are effectively assuming that GW230529 occurred at the 
end of O3. Without accounting for the extra time-volume provided by the first two weeks of O4, this introduces a 
bias in our inferred estimates of the merger rate in the mass gap (Figure 5), rate of gamma-ray bursts with NSBH 
progenitors, and total contribution of NSBH mergers to the production of heavy elements (Section 7). However this 
effect is negligible given that GW230529 occurred in the first two weeks of O4a and the detector sensitivity during this 
period was not drastically different from O3 (see Section 2), meaning that the additional time-volume unaccounted 
for in the injection sets is small. 

To generate samples from the hierarchical likelihood for our population analyses we use the DYNESTY nested sam- 
pler (Speagle 2020) as implemented in GW POPULATION (Talbot et al. 2019) for the POWER LAW + Dip + BREAK 
and NSBH-pop analyses and the Hamiltonian Monte Carlo sampler implemented in the PYMC package (Salvatier 
et al. 2016) for the BINNED GAUSSIAN PROCESS analysis. 


L EQUATION OF STATE 


To assess the effect of the EoS constraint used for source classification (Section 6) and calculation of fEM-bright (Sec- 
tion 7), we repeat these analyses using GP-EoS constraints additionally conditioned on NICER observations (Legred 
et al. 2021, 2022) of J0030--0451 (Miller et al. 2019; Riley et al. 2019; Vinciguerra et al. 2024) and J07404-6620 (Miller 
et al. 2021; Riley et al. 2021; Salmi et al. 2022, 2023). 
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We find that including the NICER observations only changes the source classification results by at most a few 
percent. This is because the constraint on the maximum neutron star mass used in the source classification analysis is 
primarily driven by the observation of massive pulsars. Changing the information included within the EoS constraint 
leads to a smaller effect than the choice of astrophysical mass and spin distribution, as discussed in Section 6. 

The inference on the fraction of NSBH systems that may have EM counterparts changes more significantly when 
NICER observations are included in the EoS constraint. This is because the GW and pulsar-only results provide 
much more support for compact neutron stars, which inhibits the probability that the neutron star is disrupted and 
suppresses the probability of an EM counterpart. The constraints including the NICER observations favor stiffer EoSs, 


enhancing the inferred EM-bright fraction, fEM-bright = 0.13 


urn and pulling the posterior peak away from zero. 
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